import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mplt
import seaborn as sns
import arviz as az
import datetime

import fit_arima

import open_data
df_temp = open_data.open()
df = df_temp[0]

ritorno = fit_arima.compute(2, 1, catene=4, samples_per_chain=2500, burnin=1000)

20:26:49 - cmdstanpy - INFO - CmdStan start processing

chain 1 |          | 00:00 Status

chain 2 |          | 00:00 Status

chain 3 |          | 00:00 Status

chain 4 |          | 00:00 Status

20:43:33 - cmdstanpy - INFO - CmdStan done processing.
20:43:33 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: code_model_namespace::log_prob: phi[1][14] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][18] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[2][23] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][4] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][6] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[1][5] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[2][7] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[1][4] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[1][2] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/code.stan', line 157, column 4 to column 34)
Consider re-running with show_console=True if the above output is unclear!

/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/fit_arima.py:41: FutureWarning: Indexing with a float is deprecated, and will raise an IndexError in pandas 2.0. You can manually convert to an integer key instead.
  y_missing_list.append({'Stazione':df.columns[v[i]],'Data':df.index[u[i]],'Samples':inference_data.posterior.w.values[:,:,i].reshape(catene*samples_per_chain)})

Tempo di computazione: 16:05

res = az.loo(ritorno['inference_data'], var_name="log_lik", pointwise=True, scale='deviance')
res

/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:803: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

             Estimate       SE
deviance_loo -14500.04   310.21
p_loo          257.53        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)     17102   99.5%
 (0.5, 0.7]   (ok)          76    0.4%
   (0.7, 1]   (bad)         13    0.1%
   (1, Inf)   (very bad)     0    0.0%

aux_plt = az.plot_khat(res, show_bins=True, figsize=(7,5))

res = az.waic(ritorno['inference_data'], var_name="log_lik", scale='deviance')
res

/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:1645: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(

Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

              Estimate       SE
deviance_waic -14522.44   309.98
p_waic          246.33        -

There has been a warning during the calculation. Please check the results.

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['phi','theta','hyper_gamma_phi_m','hyper_gamma_phi_s'], divergences=True, figsize=(15,15))

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['mh_sigma','sh_sigma','sigma'], divergences=True, figsize=(15,10))

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['y_0','betas','c'], divergences=True, figsize=(15,10))

df_risultati = az.summary(ritorno['inference_data'], var_names=['betas'], hdi_prob=0.95)
df_risultati.loc[:,'Name'] = df_temp[1].columns
df_risultati.set_index('Name', inplace=True)
df_risultati

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Name
Altitude	-0.094	0.011	-0.116	-0.073	0.0	0.0	11468.0	7201.0	1.0
Type_Industrial	0.023	0.021	-0.017	0.064	0.0	0.0	13914.0	7421.0	1.0
Type_Traffic	0.049	0.016	0.019	0.083	0.0	0.0	10027.0	7656.0	1.0
Type_Rural	-0.037	0.019	-0.074	-0.002	0.0	0.0	10059.0	7909.0	1.0
Type_SubUrban	0.000	0.016	-0.032	0.032	0.0	0.0	9407.0	7999.0	1.0
Zone_Agglomerato	-0.097	0.022	-0.140	-0.051	0.0	0.0	14890.0	8242.0	1.0
Zone_Appennino	-0.062	0.038	-0.140	0.008	0.0	0.0	11630.0	7796.0	1.0
Zone_Est	-0.088	0.013	-0.114	-0.063	0.0	0.0	11480.0	8199.0	1.0

import time
from matplotlib.patches import Rectangle

plt.figure(figsize=(10,5))
ax = plt.subplot(1,1,1)
plt.xlim(left = -0.3, right = 0.3)
pos_rows = np.arange(0,len(df_risultati.index))
plt.ylim(bottom=-0.75, top=pos_rows[-1]+1.0)
i = 1
for row in df_risultati.index:    
    lci = df_risultati.loc[row, 'hdi_2.5%']
    rci = df_risultati.loc[row, 'hdi_97.5%']
    len_ci = rci - lci
    if lci*rci <= 0:
        col = 'red'
    else:
        col = 'blue'
    ax.add_patch(Rectangle((lci, pos_rows[-i] - 0.75/2), len_ci, 0.75,
             edgecolor = 'grey',
             facecolor = col,
             fill=True,
             lw=1))
    i += 1
    
    
plt.axvline(0, 0, pos_rows[-1]+0.75, linestyle='--', color='black')
plt.grid(axis='both')
plt.yticks(ticks=np.flip(pos_rows), labels=df_risultati.index.to_list())
plt.xticks(ticks=np.linspace(-0.3,0.3,13))
plt.show()

aux_shape = ritorno['inference_data'].posterior.gamma_th.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.gamma_th.values.reshape((aux_shape[0]*aux_shape[1])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,0,:], legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,1,:], legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.c.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.c.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_mean.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_mean.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_median.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_median.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_max.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_max.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

valori_31_dic = []
for dizionario in ritorno['missing']:
    if dizionario['Data'] == '2018-12-31':
        valori_31_dic.append(dizionario['Samples'])
        
res_plot = sns.kdeplot(valori_31_dic, legend=False)

aux_shape = ritorno['inference_data'].posterior.annual_days_over_threshold.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	FEBBIO	CASTELLUCCIO	CORTE BRUGNATELLA	SAVIGNANO DI RIGO	SAN LEO	VERUCCHIO	GHERARDI	SAN PIETRO CAPOFIUME	BADIA	DELTA CERVIA	VIA CHIARINI	GIARDINI MARGHERITA	MARECCHIA	FRANCHINI-ANGELONI	PARCO BERTOZZI	VILLA FULVIA	CAORLE	SAVIGNANO	LUGAGNANO	PARCO RESISTENZA	BESENZONE	CENTO	SAN LAZZARO	PORTA SAN FELICE	S. ROCCO	SAN FRANCESCO	GAVELLO	MALCANTONE	CASTELLARANO	ROMA	DE AMICIS	FLAMINIA	PARCO MONTECUCCO	ZALAMELLA	S. LAZZARO	ISONZO	PARCO EDILCARANI	CITTADELLA	SARAGAT	GERBIDO	PARADIGNA	BOGOLESE	REMESINA	PARCO FERRARI	CENO	MONTEBELLO	TIMAVO	GIORDANI-FARNESE	GIARDINI
count	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.0000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000
mean	0.04960	0.114400	0.192100	0.412800	0.958200	12.030600	12.050700	14.477200	14.981000	15.953300	16.672700	17.848700	19.847200	20.286800	20.618500	21.483500	21.752100	22.506800	22.719000	23.330200	24.035000	24.105700	24.598200	26.081700	27.340900	29.38250	30.602700	31.168300	31.333200	32.417300	32.769600	34.475200	34.592800	36.090000	36.704000	36.984900	37.046200	37.182800	39.434300	43.215500	43.287800	44.775300	44.9319	45.022600	45.428100	53.559800	54.039400	54.699200	61.42270
std	0.24443	0.394877	0.564296	0.807503	1.381247	6.969772	7.324718	8.627695	7.891303	8.311801	9.147783	9.575478	9.657235	10.022025	10.508616	10.170687	11.005833	10.409226	10.595757	10.749322	11.231722	10.999692	11.747602	11.669209	11.580168	11.41634	12.349037	13.210716	11.478578	13.392076	13.087125	13.747657	13.948285	14.349205	13.823462	13.672453	12.717069	13.720445	14.040384	16.784083	15.845324	15.829576	15.1831	15.192524	16.608811	16.625506	17.299109	17.465378	17.56339
min	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	1.000000	2.00000	3.000000	3.000000	3.000000	3.000000	3.000000	4.000000	2.000000	0.000000	5.000000	5.000000	1.000000	3.000000	2.000000	5.000000	4.000000	4.000000	4.0000	5.000000	6.000000	8.000000	9.000000	7.000000	14.00000
25%	0.00000	0.000000	0.000000	0.000000	0.000000	7.000000	7.000000	8.000000	9.000000	10.000000	10.000000	11.000000	13.000000	13.000000	13.000000	14.000000	14.000000	15.000000	15.000000	16.000000	16.000000	16.000000	16.000000	18.000000	19.000000	21.00000	22.000000	22.000000	23.000000	23.000000	23.000000	24.000000	25.000000	26.000000	27.000000	27.000000	28.000000	27.000000	29.000000	31.000000	32.000000	33.000000	34.0000	34.000000	33.000000	42.000000	42.000000	42.000000	49.00000
50%	0.00000	0.000000	0.000000	0.000000	0.000000	11.000000	11.000000	13.000000	14.000000	15.000000	15.000000	16.000000	18.000000	19.000000	19.000000	20.000000	20.000000	21.000000	21.000000	22.000000	22.000000	23.000000	23.000000	24.000000	26.000000	28.00000	29.000000	29.000000	30.000000	31.000000	31.000000	33.000000	33.000000	34.000000	35.000000	36.000000	36.000000	36.000000	38.000000	41.000000	42.000000	43.000000	43.0000	44.000000	44.000000	52.000000	53.000000	53.000000	60.00000
75%	0.00000	0.000000	0.000000	1.000000	1.000000	16.000000	16.000000	19.000000	19.000000	20.000000	22.000000	23.000000	25.000000	26.000000	27.000000	27.000000	28.000000	29.000000	29.000000	30.000000	30.000000	31.000000	31.000000	33.000000	34.000000	36.00000	38.000000	39.000000	38.000000	40.000000	41.000000	43.000000	43.000000	45.000000	45.000000	45.000000	45.000000	45.000000	48.000000	53.000000	53.000000	54.000000	54.0000	54.000000	55.000000	64.000000	65.000000	65.000000	73.00000
max	4.00000	5.000000	12.000000	9.000000	20.000000	60.000000	66.000000	83.000000	62.000000	67.000000	73.000000	76.000000	79.000000	79.000000	89.000000	75.000000	96.000000	83.000000	88.000000	101.000000	102.000000	98.000000	106.000000	91.000000	90.000000	101.00000	97.000000	104.000000	93.000000	120.000000	95.000000	117.000000	118.000000	107.000000	119.000000	104.000000	98.000000	109.000000	107.000000	128.000000	130.000000	118.000000	160.0000	116.000000	122.000000	172.000000	154.000000	131.000000	143.00000

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_daily_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	CASTELLUCCIO	FEBBIO	CORTE BRUGNATELLA	SAVIGNANO DI RIGO	SAN LEO	VERUCCHIO	GHERARDI	BADIA	SAN PIETRO CAPOFIUME	DELTA CERVIA	VIA CHIARINI	GIARDINI MARGHERITA	MARECCHIA	FRANCHINI-ANGELONI	PARCO BERTOZZI	VILLA FULVIA	CAORLE	SAVIGNANO	LUGAGNANO	PARCO RESISTENZA	CENTO	BESENZONE	SAN LAZZARO	PORTA SAN FELICE	S. ROCCO	SAN FRANCESCO	GAVELLO	MALCANTONE	CASTELLARANO	ROMA	DE AMICIS	PARCO MONTECUCCO	FLAMINIA	ZALAMELLA	S. LAZZARO	ISONZO	CITTADELLA	PARCO EDILCARANI	SARAGAT	GERBIDO	PARADIGNA	BOGOLESE	CENO	PARCO FERRARI	REMESINA	TIMAVO	MONTEBELLO	GIORDANI-FARNESE	GIARDINI
count	10000.0	10000.0	10000.0	10000.0	10000.0	10000.000000	10000.00000	10000.0000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.00000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000
mean	0.0	0.0	0.0	0.0	0.0	0.006600	0.00990	0.0184	0.02370	0.025400	0.040100	0.051700	0.067300	0.079200	0.089600	0.09500	0.109300	0.110000	0.118100	0.131200	0.143100	0.146900	0.163900	0.195900	0.220400	0.263700	0.307300	0.328800	0.328800	0.365900	0.374100	0.424400	0.427100	0.468600	0.490900	0.500300	0.507300	0.50920	0.578100	0.645700	0.660800	0.698800	0.70590	0.71430	0.716000	0.866500	0.86950	0.872700	0.944300
std	0.0	0.0	0.0	0.0	0.0	0.080976	0.09901	0.1344	0.15212	0.157345	0.196204	0.221432	0.250553	0.270064	0.285622	0.29323	0.312031	0.312905	0.322743	0.337636	0.350192	0.354024	0.370203	0.396912	0.414537	0.440661	0.461398	0.469801	0.469801	0.481706	0.483914	0.494276	0.494682	0.499038	0.499942	0.500025	0.499972	0.49994	0.493887	0.478325	0.473461	0.458802	0.45566	0.45177	0.450959	0.340131	0.33687	0.333325	0.229353
min	0.0	0.0	0.0	0.0	0.0	0.000000	0.00000	0.0000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.000000	0.00000	0.000000	0.000000
25%	0.0	0.0	0.0	0.0	0.0	0.000000	0.00000	0.0000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	1.000000	1.00000	1.000000	1.000000
50%	0.0	0.0	0.0	0.0	0.0	0.000000	0.00000	0.0000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.00000	1.00000	1.000000	1.000000	1.00000	1.000000	1.000000
75%	0.0	0.0	0.0	0.0	0.0	0.000000	0.00000	0.0000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.00000	1.00000	1.000000	1.000000	1.00000	1.000000	1.000000
max	0.0	0.0	0.0	0.0	0.0	1.000000	1.00000	1.0000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.00000	1.00000	1.000000	1.000000	1.00000	1.000000	1.000000

aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)

<AxesSubplot: >

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_annual_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	CASTELLUCCIO	BESENZONE	PORTA SAN FELICE	ZALAMELLA	VILLA FULVIA	GHERARDI	PARCO FERRARI	PARCO MONTECUCCO	CENTO	DELTA CERVIA	BOGOLESE	CENO	ISONZO	S. LAZZARO	S. ROCCO	MALCANTONE	FLAMINIA	CITTADELLA	SAN FRANCESCO	SARAGAT	PARCO EDILCARANI	CAORLE	ROMA	FEBBIO	CORTE BRUGNATELLA	SAVIGNANO DI RIGO	SAN LEO	VERUCCHIO	BADIA	GIARDINI MARGHERITA	CASTELLARANO	PARCO BERTOZZI	PARCO RESISTENZA	DE AMICIS	VIA CHIARINI	FRANCHINI-ANGELONI	SAVIGNANO	SAN PIETRO CAPOFIUME	SAN LAZZARO	LUGAGNANO	GAVELLO	MARECCHIA	REMESINA	GERBIDO	PARADIGNA	MONTEBELLO	GIORDANI-FARNESE	GIARDINI	TIMAVO
count	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0000	10000.0000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000
mean	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0001	0.0001	0.000300	0.000400	0.000500	0.000500	0.000600
std	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0100	0.0100	0.017319	0.019997	0.022356	0.022356	0.024489
min	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000
50%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000
75%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000
max	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0000	1.0000	1.000000	1.000000	1.000000	1.000000	1.000000

aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)

<AxesSubplot: >

y_post_pred = ritorno['inference_data'].posterior.y_post_pred.to_numpy()
y_post_pred = y_post_pred.reshape(y_post_pred.shape[0]*y_post_pred.shape[1], y_post_pred.shape[2], y_post_pred.shape[3])
y_post_pred.shape

(10000, 365, 49)

import ipywidgets as widgets
from ipywidgets import HBox, VBox
from IPython.display import display

@widgets.interact(stazione = df.columns)
def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()

interactive(children=(Dropdown(description='stazione', options=('CASTELLUCCIO', 'FEBBIO', 'CORTE BRUGNATELLA',…

def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()

for stazione in df.columns:
    f(stazione)